function [W, evals] = slwhiten(type, A, varargin)
%SLWHITEN Computes the transform to whiten a sample space
%
% [ Syntax ]
%   - W = slwhiten('cov', C, ...)
%   - W = slwhiten('samples', X, ...)
%   - [W, evals] = slwhiten(...)
%
% [ Arguments ]
%   - C:        the covariance matrix (d x d)
%   - X:        the sample matrix (d x n)
%   - W:        the transform matrix
%
% [ Description ]
%   - W = slwhiten('cov', C, ...) computes the whiten transform W, such 
%     that W^T * C * W = I.
%
%   - W = slwhiten('samples', X, ...) computes the whiten transform W, such
%     that the covariance of W^T * X is I.
%
%     You can specify the following options for the computation
%     \{:
%        - preserve:       the way to determine the dimension preseved in 
%                          the transformed space.
%
%                          It should be in form of a cell array as {cri}
%                          or {cri, v}, where cri is the criteria used for
%                          decision, and v is the extra parameter.
%
%                          The following criteria is available
%                          \{:
%                             - rank:   preserve all eigenvalues not
%                                       smaller than eps(max(eigvals))
%                             - num:    preserve all the v leading
%                                       eigenvalues
%                             - ratio:  preserve all eigenvalues not
%                                       smaller than (v * max(eigvals)
%                             - energy: preserve energy of ratio not less
%                                       than v.
%                          \:}
%                          By default, the option is set to {'rank'}.
%
%        - maxrank:        the maximum possible rank of the transform.
%                          default = [], which means to determine it
%                          based on input.
%
%        - weights:        the weights of samples. By default, it is [],
%                          which means that all samples share the same
%                          weights.
%
%                          It only takes effect when computing the
%                          transform from samples.
%
%        - vmean:          The pre-computed mean vector of the samples. 
%                          default = [], which means the mean vector has
%                          not been computed.
%
%                          It can also be given as a column vector, or a
%                          zero scalar indicating zero mean vector.
%
%                          It only takes effect when computing the
%                          transform from samples.
%
%        - method:         The method to compute the transform
%                          \{:
%                             - auto:   automatically select the best
%                                       method
%                             - std:    standard way, first compute the
%                                       covariance.
%                             - svd:    use SVD to decompose the sample
%                                       matrix.
%                             - trans:  use a transposed way, that is to
%                                       solve the eigenvalues of X^T * X
%                                       instead, and then transform back
%                                       to the original domain.
%                          \:}
%                          By default, the method is 'auto'.
%     \:}
%
%   - [W, evals] = slwhiten(...) additionally returns the row vector with
%     eigenvalues associated with the columns in W.
%
% [ Remarks ]
%   - The computation of whitening transform is quite similar to the
%     training of PCA. Hence, you may see slwhiten and slpca/slcovpca 
%     share some options.
%
% [ History ]
%   - Created by Dahua Lin, on Jul 23, 2007
%

%% parse and verify input arguments

error(nargchk(2, inf, nargin));

assert(ischar(type) && (strcmpi(type, 'cov') || strcmpi(type, 'samples')), ...
    'sltoolbox:slwhiten:invalidarg', ...
    'The first argument should be a string equaling either ''cov'' or ''samples''.');
use_samples = strcmpi(type, 'samples');

if use_samples
    X = A;
    assert(isnumeric(X) && ndims(X) == 2, 'sltoolbox:slwhiten:invalidarg', ...
        'The sample matrix X should be a 2D numeric matrix.');
else
    C = A;
    assert(isfloat(C) && ndims(C) == 2 && size(C,1) == size(C,2), ...
        'sltoolbox:slwhiten:invalidarg', ...
        'The covariance matrix should be a numeric square matrix.');
end

% options

opts = struct( ...
    'preserve', {{'rank'}}, ...
    'maxrank', [], ...
    'weights', [], ...
    'vmean', [], ...
    'method', 'auto');

if ~isempty(varargin);
    opts = setopts(opts, varargin);
end

preserve = opts.preserve;
assert(iscell(preserve) && (numel(preserve) == 1 || numel(preserve) == 2), ...
    'sltoolbox:slwhiten:invalidopt', ...
    'The option preserve should be a cell array with 1 or 2 elements.');
fh_seleig = geteigpreserve(opts.preserve{:});

maxrank = opts.maxrank;
if ~isempty(maxrank)
    assert(isnumeric(maxrank) && isscalar(maxrank) && maxrank > 0 && maxrank == fix(maxrank), ...
        'sltoolbox:slwhiten:invalidopt', ...
        'The option maxrank should be a positive integer scalar.');
end

if use_samples
    
    weights = opts.weights;
    if ~isempty(weights)
        assert(isnumeric(weights) && isvector(weights) && length(weights) == size(X,2), ...
            'sltoolbox:slwhiten:invalidopt', ...
            'The option weights should be a numeric vector of length n.');
    end
    
    vmean = opts.vmean;
    if ~isempty(vmean) && ~isequal(vmean, 0)
        assert(isnumeric(vmean) && isequal(size(vmean), [size(X,1)  1]), ...
            'sltoolbox:slwhiten:invalidopt', ...
            'The option vmean should be a d x 1 vector or a zero scalar.');
    end
    
    method = opts.method;
    assert(ischar(method), 'sltoolbox:slwhiten:invalidopt', ...
        'The option method should be a string indicating its name.');        
end


%% main

% solve the eigen problem

if use_samples
    % centralize
    if isempty(vmean)
        vmean = slmean(X, weights);
    end
    if ~all(vmean == 0)
        X = bsxfun(@minus, X, vmean);
    end
    
    % normalize weights
    if ~isempty(weights)
        weights = weights / sum(weights);
    end
    
    % solve
    switch method
        case 'auto'
            [P, evals] = solve_auto(X, weights);
        case 'std'
            [P, evals] = solve_std(X, weights);
        case 'svd'
            [P, evals] = solve_svd(X, weights);
        case 'trans'
            [P, evals] = solve_trans(X, weights);
    end
    
    rank_ub = min(size(X,1), size(X,2)-1);
else
    [evals, P] = slsymeig(C);
    
    rank_ub = size(C, 1);
end

% confine with maxrank

if isempty(maxrank)
    maxrank = rank_ub;
else
    maxrank = min(rank_ub, maxrank);
end

if size(P, 2) > maxrank
    P = P(:, 1:maxrank);
    evals = evals(1:maxrank);
    evals(evals < 0) = 0;
end

% preserve

if ~isempty(fh_seleig)
    inds = fh_seleig(evals);    
    P = P(:, inds);
    evals = evals(inds);    
end

% make the final transform

evals = evals(:)';
W = bsxfun(@times, P, 1 ./ sqrt(evals));


%% The core function to solve whitening transform


function [P, evals] = solve_auto(X, w)

[d, n] = size(X);
if d <= n
    [P, evals] = solve_std(X, w);
else
    [P, evals] = solve_trans(X, w);
end


function [P, evals] = solve_std(X, w)

if isempty(w)
    C = X * X' / size(X, 2);
else
    C = bsxfun(@times, X, w) * X';
end 
[evals, P] = slsymeig(C);


function [P, evals] = solve_svd(X, w)

if isempty(w)
    [P, D] = svd(X, 0);
    evals = diag(D) .^ 2 / size(X, 2);
else
    [P, D] = svd(bsxfun(@times, X, sqrt(w)), 0);
    evals = diag(D) .^ 2;
end    


function [P, evals] = solve_trans(X, w)

if isempty(w)
    Ct = X' * X / size(X, 2);
else
    Z = bsxfun(@times, X, sqrt(w));
    Ct = Z' * Z;
end
[evals, P] = slsymeig(Ct);

if isempty(w)
    P = slnrmvecs(X * P);
else
    P = slnrmvecs(Z * P);
end

